The present invention relates to a method, apparatus, and computer program for processing digital images to reduce noise.
Many image processing noise reduction algorithms can be classified as non-linear spatial filters. Often these algorithms involve using the pixel values in a small local neighborhood surrounding the pixel of interest combined with some form of non-linear weighting and/or statistical conditions applied to the pixels in the neighborhood to derive a noise free estimate of the pixel of interest. The small local neighborhood is usually centered on the pixel of interest. For this class of noise reduction algorithms the filter size is fixed, meaning that all image pixels are processed with the same size local neighborhood. The most common shape to the local neighborhood is a rectangular region centered about the pixel of interest. Such a region can be characterized by a width and height. Usually the width and height dimensions are chosen to be symmetric.
An example of a fixed size rectangular region noise reduction algorithm is the Sigma Filter, described by Jong-Son Lee in the journal article xe2x80x9cDigital Image Smoothing and the Sigma Filterxe2x80x9d, Computer Vision, Graphics, and Image Processing, Vol. 24, 1983, pp. 255-269. This is a noise reduction filter that uses a non-linear pixel averaging technique sampled from a rectangular window about the center pixel. Pixels in the local neighborhood are either included or excluded from the numerical average on the basis of the difference between the pixel and the center pixel. Mathematically, the Sigma Filter can be represented as
qmn=xcexa3ijaijpij/xcexa3ijaijxe2x80x83xe2x80x83(1)
where:
xe2x80x83aij=1 if|pijxe2x88x92pmn| less than =xcex5,
aij=0 if|pijxe2x88x92pmn| greater than xcex5.
where pij represents the pixels in the local surround about the center pixel pmn, qmn represents the noise cleaned pixel, and xcex5 represents a numerical constant usually set to two times the expected noise standard deviation. The local pixels are sampled from a rectangular region centered about the pixel of interest.
The Sigma Filter was designed for image processing applications for which the dominant noise source is Gaussian additive noise. Signal dependent noise sources can easily be incorporated by making the e parameter a function of the signal strength. However, for both signal independent and signal dependent noise cases the expected noise standard deviation must be known to obtain optimal results. The Sigma Filter performs well on highly structured areas due to the fact that most of the image pixels in the local neighborhood are excluded from the averaging process. This leaves high signal strength regions nearly unaltered. The filter also works well in large uniform areas devoid of image signal structure due to the fact that most of the local pixels are included in the averaging process. For these regions, the Sigma Filter behaves as a low pass spatial filter with a rectangular shape. This low-pass spatial filter shape does not filter very low spatial frequency components of the noise. The resulting noise reduced images can have a blotchy or mottled appearance in otherwise large uniform areas.
Regions in images characterized by low amplitude signal modulation, or low signal strength, are not served well by the Sigma Filter. For these regions, most of the local pixel values are included in the averaging process thus resulting in a loss of signal modulation. Setting the threshold of the filter to a lower value does reduce the loss of signal, however, the noise is left mostly the same.
Another example of a fixed size non-linear noise filter was reported by Arce and McLoughlin in the journal article xe2x80x9cTheoretical Analysis of the Max/Median Filterxe2x80x9d, IEEE Transactions Acoustics, Speech and Signal Processing, ASSP-35, No. 1, January 1987, pp. 60-69, they named the Max/Median Filter. This filter separated the local surround region into four overlapping regionsxe2x80x94horizontal, vertical, and two diagonal pixels with each region containing the center pixel. A pixel estimate was calculated for each region separately by applying and taking the statistical median pixel value sampled from the regions"" pixel values. Of these four pixel estimates, the maximum valued estimate was chosen as the noise cleaned pixel. Mathematically the Max/Median Filter can be represented as                                           q            ij                    =                      maximum            ⁢                          xe2x80x83                        ⁢            of            ⁢                          xe2x80x83                        ⁢                          {                                                Z                  1                                ,                                  Z                  2                                ,                                  Z                  3                                ,                                  Z                  4                                            }                                      ⁢                  
                ⁢                              Z            1                    =                      median            ⁢                          xe2x80x83                        ⁢            of            ⁢                          xe2x80x83                        ⁢                          {                                                p                                      i                    ,                                          j                      -                      w                                                                      ,                                  …                  ⁢                                      xe2x80x83                                    ⁢                                      p                    ij                                                  ,                …                ⁢                                  xe2x80x83                                ,                                  p                                      ij                    +                    w                                                              }                                      ⁢                  
                ⁢                              Z            2                    =                      median            ⁢                          xe2x80x83                        ⁢            of            ⁢                          xe2x80x83                        ⁢                          {                                                p                                                            i                      -                      w                                        ,                    j                                                  ,                                  …                  ⁢                                      xe2x80x83                                    ⁢                                      p                    ij                                                  ,                …                ⁢                                  xe2x80x83                                ,                                  p                                                            i                      +                      w                                        ,                    j                                                              }                                      ⁢                  
                ⁢                              Z            3                    =                      median            ⁢                          xe2x80x83                        ⁢            of            ⁢                          xe2x80x83                        ⁢                          {                                                p                                                            i                      +                      w                                        ,                                          j                      -                      w                                                                      ,                                  …                  ⁢                                      xe2x80x83                                    ⁢                                      p                    ij                                                  ,                …                ⁢                                  xe2x80x83                                ,                                  p                                                            i                      -                      w                                        ,                                          j                      +                      w                                                                                  }                                      ⁢                  
                ⁢                              Z            4                    =                      median            ⁢                          xe2x80x83                        ⁢            of            ⁢                          xe2x80x83                        ⁢                          {                                                p                                                            i                      -                      w                                        ,                                          j                      -                      w                                                                      ,                                  …                  ⁢                                      xe2x80x83                                    ⁢                                      p                    ij                                                  ,                …                ⁢                                  xe2x80x83                                ,                                  p                                                            i                      +                      w                                        ,                                          j                      +                      w                                                                                  }                                                          (        2        )            
where qij represents the noise cleaned pixel, Z1, Z2, Z3, and Z4 represent the four pixel estimates, and pij represents the local pixel values. The Max/Median Filter also reduces the noise present while preserving edges. For Gaussian additive noise, the statistical median value does not reduce the noise by as great a factor as numerical averaging. However, this filter does work well on non-Gaussian additive noise such as spurious noise.
Noise is most visible and objectionable in images containing areas with little signal structure, e.g. blue sky regions with little or no clouds. The Sigma filter can produce a blotchy, or mottled, effect when applied to image regions characterized by low signal content. This is largely due to the rectangular geometric sampling of local pixels strategy. The radial region sampling strategy employed by the Max/Median Filter produces noise reduced images with less objectionable artifacts in image regions characterized by low signal content. For images with high noise content, the artifacts produced by radial region sampling strategy have a structured appearance.
Nagao and Matsuyama described an edge preserving spatial filtering technique in their publication, xe2x80x9cEdge Preserving Smoothing,xe2x80x9d in Computer Vision, Graphics, and Image Processing, Vol. 9, pp. 394-407, 1979. Nagao formed multiple local regions about the center pixel by rotating a line segment inclusion mask pivoting about the center pixel. Imagine each spatial region of pixels corresponding to a different orientation of the hour hand of a clock. For each region, the statistical variance and numerical mean are calculated. The noise cleaned pixel value is assigned as the numerical mean of the region with the lowest statistical variance. This filter does not assume a prior knowledge of the noise magnitude. If the magnitude of the inherent image structure is greater than the noise, the Nagao filter will reduce some noise while preserving edge structure. Unfortunately, the filter suffers from two problems. The size of the statistical sampling is relatively small since only one local region effectively contributes to the pixel estimation process. The other problem with this filter results from the fact that some image structure content is always degraded due to the fact that at least one region""s numerical mean replaces the original pixel value. In addition, significant artifacts (distortions to the true image structure) can occur.
U.S. Pat. No. 5,671,264, issued Sep. 23, 1997 to Florent et al., entitled xe2x80x9cMethod for the Spatial Filtering of the Noise in a Digital Image, and Device for Carrying Out the Methodxe2x80x9d, describes a variation of the Sigma Filter and Max/Median Filter. This algorithm borrows the technique of radial spatial sampling and multiple pixel estimates from the Max/Median Filter. However, the algorithm expands the number of radial line segment to include configurations with more than four segments. The algorithm uses combinations of Sigma and Median filters to form the individual region pixel estimates. These pixel estimates derived from the N regions are then combined by numerical averaging or taking the statistical median value to form the noise cleaned pixel value. A key component of this algorithm is the randomization of one of the three essential region parameters: length, orientation, and number of regions. The randomization of the filter parameters is performed on a pixel to pixel basis thus changing the inherent characteristics with pixel location. It is claimed that the randomization feature reduces the induced structured artifacts produced by the radial region geometry sampling method. The imaging application cited in U.S. Pat. No. 5,671,264 is medial x-ray imagery. This type of imagery is typically characterized by high noise content or a low signal-to-noise ratio. The structured artifacts introduced by the noise reduction algorithm are worse for low signal-to-noise ratio images.
U.S. Pat. No. 5,594,816, issued Jan. 14, 1997 to Kaplan et al. entitled xe2x80x9cComputer Based Digital Image Noise Reduction Method Based on Over-Lapping Planar Approximationxe2x80x9d describes a method of noise reduction including the steps of: a) identifying a pixel of interest; b) identifying multiple rectangular pixel regions which include the pixel of interest; c) calculating a noise free pixel estimate for each rectangular pixel region by performing a least squares fit of the pixels within the rectangular region with respect to a planar surface; d) calculating a goodness-of-fit statistical value of the planar model for the individual rectangular regions; and e) calculating a noise reduced pixel value using a weighted sum of the noise free pixel estimates based on the goodness-of-fit statistical values. While the method disclosed by Kaplan et al. involves weighting the multiple noise free estimates, the method is overly complicated and the calculation of the weighting factors is integrally related to the method of calculating the noise free pixel estimates.
The prior art methods of combining multiple pixel estimates from multiple local regions of pixels all suffer from an inability to adapt effectively to local image statistics. In structurally flat regions of images, the techniques of using the maximum of the estimates does not work well due to the fact that only one of the pixel estimates can contribute to the noise cleaning. Methods that use the arithmetic mean or median work much better in flat regions but suffer from over cleaning image regions in spatially active regions.
There is a need therefore for an improved noise reduction algorithm that operates well in both the flat and spatially active regions of an image.
The need is met according to the present invention by providing a method of processing a digital image channel to remove noise, including the steps of: identifying a pixel of interest; identifying at least two sampled local regions of pixels which include the pixel of interest; calculating a noise free pixel estimate for each sampled local region of pixels; calculating a statistical weighting factor for each sampled local region, the calculation of the statistical weighting factor being independent of the calculation of the noise free pixel estimate; and using the noise free pixel estimates and the statistical weighting factors for calculating a noise reduced pixel value.
The present invention overcomes the deficiency of using only one of the calculated noise free pixel estimates to form the noise reduced pixel value as exhibited by the Nagao et al. and Arce et al. algorithms. The present invention also overcomes the deficiency of using the median or mean methods of combining multiple noise free pixel estimates to form the noise reduced pixel as described by Nagao et al. and Florent et al. which gives too much emphasis to some of the noise free pixel estimates. The present invention allows all the calculated noise free pixel estimates to contribute in spatially flat image regions and for the noise free pixel estimates to contribute unequally in highly structured regions. This is achieved by the use of a statistical weighting factor which ranges in value for the individual sample local regions of pixels.